C:::::::::::::::::::::::::::: CHEMO_VN2 :::::::::::::::::::::::::::::::
C---- this program sets up the Jacobian matrix for the band solver
      SUBROUTINE CHEMO_VN2(ID,   !.. IN: Unit number for output
     >                    IPR,   !.. IN: # of N2(v) species
     >                      J,   !.. IN: altitude index
     >                    ALT,   !.. IN: altitude (km)
     >                    N2N,   !.. IN: N2(v) populations
     >                    NVF,   !.. IN: N2(v) populations
     >                   VKEN)   !.. IN: O+ + N2 reaction rate factor  
      IMPLICIT NONE
      INTEGER VD,IS,ID,IPR,ILJ,J,ITV
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DOUBLE PRECISION B(22),PROD(7),LOSS(7),N2N,N2RAT(6)
      DOUBLE PRECISION NV(7,VD),NVS(7,VD),NVF(7,VD),ALT
      DOUBLE PRECISION VK(6),VKE,VKEN,NTOT
      DATA VK/1.0,1.0,38.0,85.0,220.0,270.0/
      DATA N2RAT/1.0,0.2,0.02,0.005,0.002,0.0006/

      DO IS=1,IPR
        NV(IS,J)=0.1*N2RAT(IS)*N2N
      ENDDO

      DO ITV=1,4
      WRITE(6,'(3X,A,8X,A,4X,A,4X,A)') 
     >    'ITV  IS      ZR      PROD      LOSS       NV       NVS'
     >     ,'N2      DNV        VKE       eprod'
     >     ,'eloss     Oprod      Oloss    NO,O1D     VCON      N2A'
     >     ,'e*prod'
        VKE=0.0
        NTOT=N2N
  
        DO IS=IPR,2,-1
          CALL VN2PRLO(0,IS,J,IPR,NV,PROD(IS),LOSS(IS),B)
          NVS(IS,J)=NV(IS,J)
          NV(IS,J)=PROD(IS)*NV(IS,J)/LOSS(IS)
          NTOT=NTOT+NV(IS,J)
          VKE=VKE+VK(IS)*NV(IS,J)
          WRITE(ID,'(2I5,1P,22E10.2)') ITV,IS,ALT,PROD(IS),LOSS(IS),
     >     NV(IS,J),NVF(IS,J),N2N,(NV(IS,J)-NVS(IS,J))/NVS(IS,J)
     >    ,VKE/NTOT,B(1),B(2),B(7),B(8),B(9),B(12),B(13),B(14)
c     >     ,VKE/(NV(1,J)+NV(2,J)+NV(3,J)+NV(4,J))
        ENDDO
      ENDDO
        NTOT=0.0
        VKE=0.0
        DO IS=1,IPR
          NTOT=NTOT+NV(IS,J)
          VKE=VKE+VK(IS)*NV(IS,J)
        ENDDO
        VKEN=VKE/NTOT
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE VN2PRLO(ILJ,I,JP,IPR,NV,PROD,LOSS,B)
C--------this subroutine evaluates the production and loss for the various
C--------states of vib n2. see newton et al jgr 1974 p3817.
      IMPLICIT REAL(A-H,L,N-Z)
      INTEGER VD
      PARAMETER (VD=181)    !.. VD = dimension of vertical grid
      DOUBLE PRECISION PROD,LOSS,B,STR,RCON,F,NV(7,VD),NVS,NNO,N2D,O1D,
     > ON,NR,ZR,TN,GR,BG,VDT,TE,VK,PXN2A,PEEX,PRN2P,O2N,TIZ,ZNE,HE,N2P
     > ,N4S,HN,RTS,N2N,CHIV,DH
      DOUBLE PRECISION TE2,TE3,AIJ,AJI,AIJN,AJIN,VCON,SQRTN,NU,P1001,
     > P1NU,EXPTN,C01,C01ON,TE1
      DIMENSION B(22),RTS(199)
      COMMON/STAW/STR(12,401),RCON(401)
      COMMON/NDM/NNO(VD),N2D(VD),O1D(VD),N4S(VD),N2N(VD),ON(VD),
     > CHIV(VD)
      COMMON/ND1/NR(2,VD),ZR(VD),TN(VD),GR(VD),BG(VD),VDT,TE(VD),DH,
     > JTER
      COMMON/PXS/VK(6),PXN2A(VD),PEEX(7,VD),PRN2P(7,VD),JVMAX
      COMMON/ND2/O2N(VD),TIZ(VD),ZNE(VD),HE(VD),N2P(VD),HN(VD)
      COMMON/PEE/P(4,6,6),G(6),HCOK
      DATA IFRST/0/
C------ the following data(g @ pc) are related to thermal electron excitation
C------ of vib n2. the g's are used to calculate the dexcitation.
C------ the p parameters are given by newton et al p3817 for electron exCit.
C
C----- the b's are used to sum the v-v exchange between diff states. however
C----- b(1) @ b(2) are used to sum elec excit. nu=coll freq eqtn(18). p1001=
C----- dimensionless exchange transition probability for n2 trans 1<--->0
C----- c01= tras-vib energy exchange prob. for colls between o and n2. however
C----- the results from mcneal et al jgr 1974 are used for reaction rate
      IF(IFRST.EQ.0)CALL INITP(IFRST)
      DO  5 IS=1,22
  5   B(IS)=0.0
      CALL RATS(JP,TE(JP),TIZ(JP),TN(JP),RTS)
      SQRTN=DSQRT(TN(JP))
      NU=1.482094E-11*SQRTN
      P1001=2.600322E-6*TN(JP)*EXP(91.5/TN(JP))
      P1NU=P1001*NU
      EXPTN=EXP(3380/TN(JP))
      C01=TN(JP)**2.87*2.4685E-22/EXPTN
      C01ON=C01*ON(JP)
      TE1=1.0/TE(JP)
      IF(TE(JP).LT.400) TE1=1.0/400.0
      TE2=TE1*TE1
      TE3=TE2*TE1
C
C++++++ compts of prod and loss eqtns 16 @ 17 of newton. note that indices
C++++++ in coefficients of production and loss terms
C+++++++ appear reduced by 1, because the fortran indices start from 1, not 0
      DO 40 J=1,IPR
      IF(I.EQ.J) GO TO 50
      IF(P(1,I,J).NE.0.0)
     >  AIJ=EXP(P(1,I,J)+P(2,I,J)*TE1+P(3,I,J)*TE2+P(4,I,J)*TE3)
      IF(P(1,J,I).NE.0.0)
     >  AJI=EXP(P(1,J,I)+P(2,J,I)*TE1+P(3,J,I)*TE2+P(4,J,I)*TE3)
      AIJN=AIJ
      AJIN=AJI
      IF(P(1,I,J).EQ.0) AIJN=AJI*EXP(+(G(J)-G(I))*HCOK*TE1)
      IF(P(1,J,I).EQ.0) AJIN=AIJ*EXP(+(G(I)-G(J))*HCOK*TE1)
C----- zne(jp) contains ne from sub prog prln
      IF(ILJ.NE.4) B(1)=B(1)+AIJN*NV(J,JP)*(ZNE(JP))
      IF(ILJ.NE.4) B(2)=B(2)+AJIN*(ZNE(JP))*NV(I,JP)
      IF(ILJ.EQ.4.AND.J.LT.I) B(1)=B(1)+AIJN*NV(J,JP)*(ZNE(JP))
      IF(ILJ.EQ.4.AND.J.LT.I) B(2)=B(2)+AJIN*(ZNE(JP))*NV(I,JP)
 50   IF(I.EQ.1) GO TO 60
      IF(J+1.GT.IPR) GO TO 60
C------ if i-j=1 quanta are exchanged but no net change in any density
      IF(I-J.EQ.1) GO TO 60
      B(5)=B(5)+P1NU*(I-1)*J*NV(I-1,JP)*NV(J+1,JP)
      B(3)=B(3)+P1NU*(I-1)*J*NV(J,JP)*NV(I,JP)
 60   IF(J.EQ.1) GO TO 40
      IF(I+1.GT.IPR) GO TO 40
C------ if j-i=1 quanta are exchanged but no net change in any density
      IF(J-I.EQ.1) GO TO 40
      B(4)=B(4)+P1NU*(J-1)*I*NV(J,JP)*NV(I,JP)
      B(6)=B(6)+P1NU*(J-1)*I*NV(I+1,JP)*NV(J-1,JP)
 40   CONTINUE
C
C
C******* b(7)=prod due to t-v between o and n2
      IF(I.NE.1) B(7)=C01ON*((I-1)*NV(I-1,JP)+I*EXPTN*NV(I+1,JP))
      IF(I.EQ.1) B(7)=C01ON*I*EXPTN*NV(I+1,JP)
C******* b(8)=loss due to quenching of n2 by o
      IF(I.LE.9) B(8)=C01ON*((I-1)*EXPTN+I)*NV(I,JP)
C******* b(9)=prod thru no+n-->n2*+o & o(1d)+n2--->. & n2(a)
C------ for pno @ po1d see waite et al pss 1979 p901.
      IF(I.EQ.5) B(9)=RTS(9)*N4S(JP)*NNO(JP)
      IF(I.EQ.3) B(9)=RTS(33)*O1D(JP)*NV(1,JP)
C------ production due to the electronic excitation of n2(a)
      IF(I.EQ.1) B(15)=1.3*PXN2A(JP)*NV(1,JP)
      IF(I.EQ.IPR) B(13)=1.3*PXN2A(JP)*NV(1,JP)
C------ quenching of n2* by o2 .......
C...      if(i.ne.1) b(10)=1.0e-13*o2n(jp)*nv(i,jp)
C...      b(11)=1.0e-13*o2n(jp)*nv(i+1,jp)
C------ loss by o+ + n2* --> o + n2+
      JCON=J
      IF(MINJ.GT.1) JCON=MAXJ+1-J
      VCON=STR(7,JCON)
      B(12)=VCON*RTS(3)*NV(I,JP)*NR(1,JP)
C------ direct excitation by photoelectrons
      IF(I.NE.1) B(14)=PEEX(I,JP)*NV(1,JP)
      IF(I.EQ.1) B(16)=PEEX(I,JP)*NV(1,JP)

      B(1)=B(1)/1.75
c      B(2)=1.75*B(2)
      
C------ sum the cmpts of prod and loss
      PROD=B(1)+B(5)+B(6)+B(7)+B(9)+B(11)+B(13)+B(14)
      LOSS=B(2)+B(3)+B(4)+B(8)+B(10)+B(12)+B(15)+B(16)
      IF(I.EQ.1) PROD=0.0
      IF(I.EQ.1) LOSS=0.0
      RETURN
      END
